library(ggplot2)
half_life <- 1.387 # half life of 10Be in Million years
initial_ratio <- 13.87e-8 # Initial Authigenic measured in Gorongosa
AMS_limit <- 1e-15 # technical limit of AMS

decay_rate <- function(t){
  return(exp(-log(2) * t / half_life))
}

times <- seq(0,50,.1) # time range in Million years
decay <- decay_rate(times) 

data <- data.frame(time = times, decay = decay)

x <- 19.5
y_interp <- approx(data$time, data$decay, x)$y
point_data <- data.frame(time = x, decay = y_interp)

y <- AMS_limit / initial_ratio
x_interp <- approx(data$decay, data$time, y)$y
vline_data <- data.frame(time = x_interp, decay = y)

ggplot(data, aes(x = time, y = decay)) + 
  geom_line(size = 2) +
  geom_point(data = point_data, color = "blue", size = 4) +
  geom_text(data = point_data, aes(label = round(decay,10)), vjust = -1, color = "blue") +
  ggtitle("Age limit for 13.87e-8 Aut. Init. ~ 37.5 Ma") +
  xlab("Time (Ma)") +
  ylab("Fraction remaining") +
  scale_x_continuous(breaks = seq(0,50,2)) +
  scale_y_continuous(label = scales::percent) +
  theme_minimal() + theme(text = element_text(size = 18)) +
  geom_vline(xintercept = vline_data$time, linetype = "dashed",
             colour = "red", size = 2) 

ggplot(data, aes(x = time, y = decay)) + 
  geom_line(size = 2) +
  geom_point(data = point_data, color = "blue", size = 4) +
  geom_text(data = point_data, aes(label = round(decay,10)), vjust = -2, color = "blue") +
  geom_text(data = point_data, aes(label = time), vjust = 2, color = "black") +
  ggtitle("Plot in log, red dash age is where AMS detection limit 1e-15 intercepts") +
  xlab("Time (Ma)") +
  ylab("Fraction remaining (log scale)") +
  scale_x_continuous(breaks = seq(0,50,2)) +
  scale_y_continuous(trans = "log", label = scales::percent) +
  theme_minimal() + theme(text = element_text(size = 18)) +
  geom_vline(xintercept = vline_data$time, linetype = "dashed",
             colour = "red", size = 2) 

Authigenic plots

# Initial Authigenic values seem to range from -6 to -12 in the literature
initial <- 10^(seq(-4,-15,length.out = length(times)))
data$initial <- initial
detected <- expand.grid(initial = data$initial, decay = data$decay)
detected$detected <- detected$initial * detected$decay

# Set threshold and range for z values
threshold <- 1e-15 # AMS limit is 1 ppq
range <- 0.1e-15
detected$time <- rep(times, each = length(decay))

limitset <- detected[detected$detected > 1e-15,]


library(plotly)
# Create a scatter plot with x, y, and z variables
p1 <- plot_ly(limitset, x = ~initial, y = ~decay, z = ~time, type = 'scatter3d',
             mode = 'markers', color = ~log(detected))  %>%
  layout(scene = list(xaxis = list(title = 'Initial Authigenic', type = 'linear', tickformat = ".2e"),
                      yaxis = list(title = 'Decay ratio (%)', type = 'linear'),
                      zaxis = list(title = 'Time (Million years)', type = 'linear')))

# Show the plot
p1
# Create a scatter plot with x, y, and z variables
p2 <- plot_ly(limitset, x = ~initial, y = ~decay, z = ~time, type = 'scatter3d',
             mode = 'markers', color = ~log(detected))  %>%
  layout(scene = list(xaxis = list(title = 'Initial Authigenic (log)', type = 'log', tickformat = ".2e"),
                      yaxis = list(title = 'Decay ratio (%)', type = 'linear'),
                      zaxis = list(title = 'Time (Million years)', type = 'linear')))

# Show the plot
p2